home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_100 / 167_01 / spline.doc < prev   
Text File  |  1985-11-13  |  21KB  |  534 lines
  1. /* 
  2. HEADER:     CUG
  3. TITLE:        Cubic Spline Functions Theory;
  4. VERSION:    3.00;
  5. DATE:        09/26/86;
  6. DESCRIPTION:    "Mathematical background and development of
  7.         equations used in SPLINE, a full emulation of
  8.         Unix's 'spline' utility.";
  9. KEYWORDS:    spline, graphics, plot, filter, UNIX;
  10. SYSTEM:
  11. FILENAME:    SPLINE.DOC;
  12. WARNINGS:
  13. CRC:        xxxx
  14. SEE-ALSO:    SPLINE.C;
  15. AUTHORS:    Ian Ashdown - byHeart Software;
  16. COMPILERS:
  17. REFERENCES:    AUTHORS: R.W. Hamming;
  18.         TITLE:     Numerical Methods for Scientists and
  19.              Engineers, 2nd Edition
  20.              pp. 349 - 355, McGraw-Hill (1973);
  21.         AUTHORS: Forsythe, Malcolm & Moler;
  22.         TITLE:     Computer Methods for Mathematical
  23.              Computations
  24.              pp. 70 - 83, Prentice-Hall;
  25. ENDREF
  26. */
  27.  
  28.  
  29. byHeart Software
  30. 620 Ballantree Road
  31. West Vancouver, B.C.
  32. Canada V7S 1W3
  33. September 13th, 1986
  34.  
  35.  
  36.         Curve Fitting With Cubic Splines
  37.         --------------------------------
  38.  
  39.              Ian E. Ashdown
  40.             byHeart Software
  41.  
  42.  
  43. ========================================================================
  44.  
  45. Abstract: CURVE FITTING WITH CUBIC SPLINES
  46. ------------------------------------------
  47.  
  48.     A rigorous development of the mathematical theory of cubic
  49. splines and their use in curve fitting in two and three dimensions.
  50. Emphasis is on the implementation of efficent computer algorithms to
  51. calculate coefficients of both nonperiodic and periodic splines using
  52. Gaussian Elimination and sparse matrix representations. Accompanying
  53. program listing presents C source code for a full emulation of the UNIX
  54. (registered trademark Bell Laboratories) SPLINE curve interpolation
  55. utlity.
  56.  
  57. ========================================================================
  58.  
  59.  
  60.  
  61.     Before computer-aided drafting workstations completely replace the
  62. draftsman's pencil and paper, let's examine one of his tools: the
  63. spline. Presented with data in the form of points on an x-y plane, the
  64. draftsman will use a spline - a flexible strip of metal or plastic - to
  65. draw a smooth curve between them.
  66.  
  67.     The technique is very simple. After plotting the data on a sheet of
  68. paper, an appropriately sized spline is held in place at these points
  69. (referred to as "knots") with weights or pins. The draftsman then traces
  70. the curve formed by the spline. For any given set of knots, the curve
  71. generated is independent of the spline chosen, and is thus exactly
  72. reproducible.
  73.  
  74.     From mechanical engineering, elementary beam theory shows that if
  75. the spline is not too severely stressed, it will conform to a curve
  76. described by a set of cubic polynomial equations, one between each pair
  77. of adjacent knots. Adjacent polynomials meet at their common endpoints
  78. (the knots), and their slopes and rates of curvature at these points are
  79. equal. Stated in mathematical terms, these polynomials join continuously
  80. at the knots with continuous first and second derivatives.
  81.  
  82.     Knowing this, we can develop a mathematical model of the draftsman's
  83. splines, and from this model construct a computer program for
  84. interpolating a smooth curve between a set of knots. With a bit of
  85. care in choosing algorithms, such a program can quickly and accurately
  86. generate a curve for a thousand or more data points on the smallest of
  87. personal computers. It can even be adapted to interpolate a smooth
  88. surface between points plotted in three dimensions.
  89.  
  90.     Developing the model involves basic calculus and matrix theory. If
  91. you are unfamiliar with such mathematics, rest assured that the
  92. resultant algorithms are very easy to program, and using a cubic spline
  93. program requires no understanding of the underlying mathematical theory.
  94. Give the program a set of knots and it will dutifully interpolate a
  95. smooth curve in all (well, almost) cases.
  96.  
  97.     Why then discuss the mathematics of cubic splines at all? There are
  98. two answers. One is that seeing how the algorithms are developed gives
  99. you the confidence to use them. The other is that there may be cases
  100. where the algorithms will not perform exactly as desired. Knowing their
  101. theory may enable you to create a modified algorithm to fit the problem
  102. at hand.
  103.  
  104. A Simple Explanation
  105. --------------------
  106.  
  107.     Whether you understand calculus and matrices or not, it helps to
  108. have in mind a mental image of what you are trying to accomplish. Even
  109. if you aren't comfortable with the mathematics, a conceptual model will
  110. serve to illustrate the principles involved.
  111.  
  112.     A cubic polynomial equation is nothing more than a simple algebraic
  113. equation of the form:
  114.  
  115.               3     2
  116.         y = ax  + bx  + cx + d
  117.  
  118. where a,b,c and d are constants. If you plot the knots on a grid and
  119. actually bend a draftsman's spline to fit between them, each section of
  120. the spline between adjacent knots can be matched by a curve generated
  121. by the above equation ... if the appropriate constants are chosen. The
  122. trick is to find those constants for each section of the spline.
  123.  
  124.     Common sense says three things about these curves. One is that they
  125. should join at the knots. Next, the slopes of the curves should be the
  126. same where they join. Third, the curvature of the curves where the join
  127. should be the same. The first is obvious - the curves must form a
  128. continuous line. The second and third are more intuitive, but a few
  129. moments thought should convince you. A semirigid object such as a
  130. spline does not permit abrupt changes in angle or curvature as it is
  131. bent around a rigid pin.
  132.  
  133.     Stating these ideas mathematically will enable you to calculate the
  134. constants for all the cubic polynomial equations needed to model the
  135. spline.
  136.  
  137. A Rigorous Explanation
  138. ----------------------
  139.  
  140.     Beginning in this section, the mathematical theory of cubic splines
  141. will be rigorously developed. If this approach does not appeal to you,
  142. simply skim the text for the information on how to implement the spline
  143. algorithms, then examine the accompanying program listing. The program
  144. header provides all the information you will need to use the program.
  145.  
  146.     Starting with a set of data points (the knots) stated as horizontal
  147. coordinates x[i] (i = 1,...,n) and corresponding vertical coordinates
  148. y[i], the curve-to-be is defined as the composite function:
  149.  
  150.     y(x) = f[i](x) for x[i] <= x < x[i+1], i = 1,...,n-2
  151.                    and x[n-1] <= x <= x[n]
  152.  
  153. where each function f[i](x) is a cubic polynomial as described above.
  154. In other words, y(x) is really a set of functions, each of which is
  155. defined over an interval between two adjacent knots at (x[i],y[i]) and
  156. (x[i+1],y[i+1]).
  157.  
  158.     Furthermore, let's define y'[i] and y"[i] as the first and second
  159. derivatives of y(x) at x = x[i]. Knowing that the set of functions
  160. f[i](x) must join at their endpoints (the knots of the spline) and also
  161. that their first and second derivatives are continuous at these points,
  162. we have the following continuity conditions:
  163.  
  164.     f[i](x[i]) = y[i]        i = 1,...,n-1
  165.     f[i-1](x[i]) = y[i]        i = 2,...,n
  166.     f'[i-1](x[i]) = f'[i](x[i])    i = 2,...,n-1
  167.     f"[i-1](x[i]) = f"[i](x[i])    i = 2,...,n-1
  168.  
  169.     Since each function f[i](x) is a cubic polynomial, it follows that
  170. its second derivative is a linear function (a straight line) between its
  171. endpoints. If we define:
  172.  
  173.     h[i] = x[i+1] - x[i]
  174.  
  175. then linear interpolation gives us:
  176.  
  177.            y"[i] * (x[i+1] - x) + y"[i+1] * (x - x[i])
  178.     f"[i](x) = -------------------------------------------
  179.                      h[i]
  180.  
  181.     Integrating this equation twice and selecting the constants of
  182. integration such that the continuity conditions are satisfied, we can
  183. derive the following interpolation equation:
  184.  
  185.           y[i] * (x[i+1] - x) + y[i+1] * (x - x[i])
  186.     f[i](x) = ----------------------------------------- -
  187.                    h[i]
  188.  
  189.               2  +-       +-
  190.           h[i]     |       | (x[i+1] - x)
  191.           ---- * | y"[i] * |------------- -
  192.            6     |       |    h[i]
  193.              +-       +-
  194.  
  195.           +-          -+ 3 -+           +-
  196.           | x[i+1] - x |    |           | (x - x[i])
  197.           | ---------- |    | +  y"[i+1] * | ---------- -
  198.           |    h[i]    |    |           |   h[i]
  199.           +-          -+   -+           +-
  200.  
  201.           +-        -+ 3 -+  -+
  202.           | x - x[i] |      |   |
  203.           | -------- |      |   |
  204.           |   h[i]   |      |   |
  205.           +-        -+     -+  -+
  206.  
  207.     for i = 1,...,n-1
  208.  
  209.              Interpo